Temporal and structural characteristics of a two dimensional gas of hard needles 
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We have simulated the dynamics ol a 2D gas ol hard needles by event-oriented molecular dy- 
namics. Various quantities namely translational and rotational diffusion constants and intermediate 
self scattering function have been explored and their dependence on density is obtained. Despite 
absence of positional ordering, the rotational degree of freedom behaves nontrivially. Slowing down 
is observed in the angular part of the motion. It is shown that above a certain density the rotational 
mean squared displacement exhibits a three stage regime including a plateau. 



I. INTRODUCTION 



Molecular dynamics (MD) simulation of anisotropic 
hard objects such as ellipsoids and spherocylinders has 
been the subject of exploration in past decades [1]. It 
has been shown the alignment of non-spherical molecules 
can lead to a diversity of phases, mainly orientational in 
nature, in liquid crystals [l|-|3|. Despite the profound in- 
sight obtained via tremendous Monte Carlo (MC) sim- 
ulations of the static phases, many dynamical, trans- 
portational and structural properties namely kinetic ar- 
rest and glassy behaviours have only been poorly under- 
stood (4|]. In spite of employment of other simulation 
techniques like Brownian dynamics MWi and theoretical 
approaches such as kinetic theory [8, <^, density func- 
tional [10I4T2I and hydrodynamics equations approach 
p^ . event-oriented MD remains as an efficient tool for 
probing the dynamical aspects of hard gases of non- 
spherical objects. Among elongated and anisotropic hard 
bodies, infinitely thin needles have received quite notable 
attention [l[ . The first MD simulation attempt were car- 
ried out for a three dimensional gas of hard needles by 
Frenkel and Maguire [H [ll| and Magda et al. [11]. In 
these investigations, various temporal auto-correlations 
were explored and compared to predictions of Enskog 
and Doi-Edward theories. The physics of hard needle gas 
is sensitive to dimensionality. By extensive MC simula- 
tions in 2D, Frenkel and Eppenga showed the existence 
of quasi long range order in the system at high densi- 
ties [17.]. It was argued that this 2D system undergoes 
a Kosterlitz-Thouless phase transition [l^. Resurgence 
of interest in the dynamical properties of hard needle gas 
was sparked by papers of Renner et al. [l9| and Obukhov 
et al. (20j who introduced a rotator model. Their simple 
ideal glass former can mimic the basic dynamical prop- 
erties of an orientational glass. This rotator model was 
shown to exhibit orientational glassy behaviour like its 
positional counterpart i.e., hard-sphere system [2l|. A 
natural question to ask is whether the orientational glassy 
behaviour survives if one releases the positional degrees 
of freedom. Chrzanowska et al. carried out the first MD 
simulation of a two dimensional hard needle gas [13, [2^ 
and mainly studied velocity auto correlations whereas 



the transport properties remained unexplored. Subse- 
quently transport coefficients like self-diffusion and shear 
viscosity of a 3D hard needle gas were studied within 
MD approach by Mukoyama et al. [13, HI] . We wish 
to note that the hard needle rotator model is amenable 
to comparison with experiments carried out via neutron 
scattering and has shown to successfully model the ex- 
perimentally observed phenomenon of orientation glassy 
dynamics f26]. Recent investigations on 2D hard needle 
gas involve deposition-evaporation dynamics [131 and re- 
striction of needles center of masses [11] . Here we focus 
on the transportation features to illuminate the interplay 
of translational and rotational degrees of freedom in 2D 
systems of highly anisotropic hard objects. 



II. 



DESCRIPTION OF THE PROBLEM 



Consider a system comprising of infinitely thin hard 
needles with mass m and length / which are restricted 
to move in the two dimensional x — y plane. The inter- 
action between needles is assumed to be hard core. We 
work in the reduced units in which m and I are taken 
unity. Upon a collision between two needles, the centre 
of mass (CM) velocity and the angular velocity around 
the z axis through the CM will change. We assume the 
collisions are elastic and frictionless therefore the energy 
(entirely kinetic) is conserved after a collision. The kinet- 
ics of such collision has been investigated in detail in [l^ . 
We study the dynamics of this 2D needle gas within the 
collision-oriented MD approach. In the collision-oriented 
MD the evolution of system takes place from collision to 
collision. The main task in MD of hard objects is to find 
the collision time between two apart objects [11]. A nu- 
merical scheme has been originally introduced in (sol . [3l| . 
The method is based on finding a geometrical condition 
for overlapping between objects. This condition is given 
as Hi,H2 < |] sin(6')] where 6 is the angle between nee- 
dles. H2 is the distance of needle 2 CM to the line of 
needle 1 and reverse. To find the next collision time we 
move all the needles in small time steps until the needles 
become very close to each other that can be regarded 
as touching. Then we update the velocity and angular 
velocity of the colliding needles and proceed to the next 
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FIG. 1: Fig.l: Collision frequency F vs p. The dependence 
is linear like 3D gas of hard needles investigated in [isl . [T^ . 
The slope changes around p = 5.5. Line is for helping eyes 




FIG. 2: Fig.2: System pressure P (reduced units) vs p. Pres- 
sure slope changes around p = 5.5. Line is for helping eyes. 



collision event. 



III. MOLECULAR DYNAMICS SIMULATION 

A. Static properties and velocity autocorrelations 

We have simulated the dynamics of a 2D gas of hard 
needles with the method explained in section II. Periodic 
boundary condition is imposed. The size of our simula- 
tion box is set to L — 7 . The number density is defined 
as p — where N is the number of needles. The total 
energy of the system E (entirely kinetic) is divided into 
two segments of translational Etrs and rotational Erot- 
We remark that due to lack of any energy scale in the 
potential energy between particles, the temperature T 
appears as an overall multiplicative factor in thermody- 
namic quantities such as pressure and free energy. Thus 
the state of the system trivially depends on temperature 
and hence energy. Nevertheless E determines the time 
scale T. 



In the thermal unit one has r = I 



have taken the energy per particle ^ 



We 



1.5 in our sim- 
ulations which gives UbT = I. A useful quantity is the 
average number of collision per unit time T that a nee- 
dle experiences. In and [11] it is argued that T has 
a linear dependence on number density p in a 3D gas. 
According to Frenkel and Maguire arguments in 3D, F 
scales as follows: F ~ p. In fact the Doi-Edward 

theory predicts such linear behaviour as well. We have 
computed the dependence of F on p in our 2D model. 
Figure (1) depicts that F has, analogous to the 3D case 
studied in jl^ J_6], a linear dependence on p. Note the 
slope change about p ~ 5.5. 

The equation of state i.e., dependence of P on the p is 
shown in figure (2). In fact P does not linearly increases 
with p. The slope changes about the same value p = 
5.5 which is in agreement with MC results of Frenkel 




FIG. 3: Fig.3: Dependence of nematic order parameter S on 
p. The transition to nematic ordered phase is expected to be 
a finite size effect. Line is for helping eyes. 



and Eppenga [17[. This marks that from translational 
viewpoint the system is totally distinct from ideal gas. 

Next we turn to order parameter. The orientational 
nematic order parameter S is defined as follows: 



1 ^ 

^=]^<E cos(20,-20,)). 



(1) 



In which 0i is the director angle of z-th needle with 
positive X axis and the average is taken over trajectories. 
Dependence of 5* on p is sketched in figure (3). 

Our simulation shows a large value for the nematic or- 
der parameter S at high densities which corresponds to a 
nematic phase. This result is expected to be a finite-size 
effect. We argue that the possibility of truly phase transi- 
tion to a nematic phase is not excluded in 2D. More con- 
cisely, the conditions of the Mermin- Wagner theorem are 
not satisfied here due to the fact that inter-molecular po- 



3 



tential between needles V{r, 9) is not separable. Despite 
Monte Carlo simulations have not shown the existence of 
such isotropic-nematic transition in a 2D gas of needles 
17 1 , a quasi long range order of Kosterlitz-Thouless type 
17, [H I was shown to persist in the system. The subject 
of temporal velocity auto correlation function has been 
the extensively studied in flF, Tg'] within the framework 
of MD and recently in [22, 23] both in event-oriented 
MD approach and analytically in the framework of En- 
skog kinetic theory. We have explored the autocorrela- 
tion between longitudinal and transverse decomposition 
of velocity. These quantities are defined as follows: 



Cu{t) 



{v{t).u{0)v{Q).u{Q)). 



(2) 




0.1 0.2 

t (reduced unit) 



(3) 



u denotes the unit vector along the needle orientation 
and the matrix P = 1 — u{0)u*{0) is the projection op- 
erator. In figure (4) we exhibit the temporal dependence 
of C||(t) and C±{t) for various densities. 

We refer the readers to [11] for a detailed discussion 
on velocity autocorrelations. We now consider the tem- 
poral auto correlation of the second order angular order 
parameter C2{t). This quantity is defined as follows: 
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C2{t) - {P2{u{t).m))- (4) 

In which P2 is the second order Legendre polynomial. 
Figure (5) shows the temporal dependence of C2{t). 

In low densities we see a fast decay which is attributed 
to the fluid like behaviour. By increasing the density the 
temporal behaviour becomes slow and two characteris- 
tics time scale emerge. This confirms the nontrivial role 
played by the angular degree of freedom. We emphasize 
that simulation on larger system is needed to ensure the 
persistence of such large correlations. The intermediate 
self scattering function Fs{q,t) is shown in figure (6). 

B. Structural properties 

We have obtained and explored various structural 
quantities. The radial distribution function g{r) (not 
shown) is featureless for r > 1 and regardless of the den- 
sity value it approaches to unity without showing any sig- 
nificant amplitude fluctuations. This demonstrates that 
the gas posses no positional ordering. However, the ex- 
istence of excluded volume effect discriminates its posi- 
tional features to ideal gas. Next we consider the angular 
spatial correlation function 52 (^)- This quantity is de- 
fined as [13: g2ir) = (cos(2je'(r) - 61(0)])). The average is 
over all needle pairs having CM to CM distance r. Figure 
(7) plots the dependence oi g2{r) vs r. 



FIG. 4: Fig.4: Temporal dependence of longitudinal and 
transverse components of velocity ACF. Longitudinal com- 
ponent exhibits a slower decay which is due to direction of 
impulsive force between needles. 




t (reduced time) 

FIG. 5: Fig. 5: Temporal dependence of C2 at various densi- 
ties (semi- log scale). 



FIG. 6: Fig. 6: Time dependence of the intermediate self 
scattering function Fs{q,t) at g = 3. 



FIG. 8: Fig. 8: Dependence of algebraic decay exponent 772 
on p. 




FIG. 7: Fig. 7: Dependence of g2{r) on r for various densities. 



For small densities, g2('') rapidly approaches zero. 
This decrease is faster than algebraic. Contrary, for high 
density p > 6 we have a slow decrease which can be in- 
dicative of slow dynamics and angular structural arrest. 
Taking the decay form algebraic as g2i'r) ~ r~^^ we have 
found out by fitting a curve the decay exponents at vari- 
ous densities. Figure (8) shows our exponents which are 
compared to those of Frenkel and Eppenga [17]. At high 
densities MD and MC results are in good agreement while 
quite notable differences are seen at low densities. 
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The term in the bracket is the translational mean- 
square displacement (MSD) ((Ar)^). Here the spatial 
dimension d is two and the average is doubly taken over 
trajectories of needles' CM and time origins. Figure (9) 
sketches the time dependence of the translational MSD 
(see journal printed version to see figure nine). After a 
ballistic regime, one recovers normal diffusive behaviour. 
This confirms that from translational viewpoint, the sys- 
tem resembles an structureless gas. Figure (10) plots a 
CM trajectory of a typical needle at two densities p = 3.9 
and p — 7.9 : 

Three distinctive kinds of motion can be identified: dif- 
fusion, channeling and entanglement among cages formed 
by neighbouring needles. At low densities the motion re- 
sembles a normal diffusion. At higher densities a chan- 
neling type of motion emerges and at higher densities the 
topological constraints in 2D leads to entanglement. Let 
us now explore the rotational diffusion Drot- In sharp 
contrast, Drot seems to be entirely of different nature. 
We have: 
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Drot — 



2CNt 



(6) 



IV. TRANSPORT PROPERTIES 

Now we report our results for the transport properties 
and orientational structure of the system. We begin with 
translational diffusion coefficient Dtrs defined as follows: 



Analogously the bracketed term denotes the angular 
mean square displacement ((A0)^) and C is the num- 
ber of angular degrees of freedom (here ^ = 1). Figure 
(11) sketches the time dependence of {{A6)^) (see journal 
printed version to see figure 11). A significant difference 
is seen compared to the translational diffusion. For p > 8 
the rotational MSD exhibits a three-stage regime which 
can be attributed to angular glassy dynamics. This is 
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FIG. 9: Fig. 10: CM trajectory of a needle at p = 3.9 during 
1586 time unit and p — 7.9 during 2329 time unit. Note 
the unequal axes scaling. The true trajectories are elongated 
along the axis having larger scale. 

consistent to the density proposed in [l3|- The possi- 
bility of angular glassy behaviour has been earlier ex- 
plored by Renner et al. [l^. It was shown above a pc, 
the angular auto correlations exhibit slow dynamics and 
multi step relaxation which can be attributed to glassy 
behaviour. Our findings is supportive of the existence 
of angular glassy behaviour even when the translational 
degree of freedom is released. We stress that simulations 
with a larger system size is crucially needed to confirm 
this conclusion. Finally we discuss the dependence of 
translational and rotational diffusion constants on p (fig- 
ures 12 and 13). 

For p < 5 translational diffusion constant Dtrs shows a 
decreasing behaviour. This is in accordance with the En- 
skog kinetic theory !15]. In fact at low densities successive 
binary collisions are uncorrelated therefore both Dtrs and 
Drot would be inversely proportional to p [1]. For large 
densities Dtrs shows an increasing trend which is in qual- 
itative agreement with the 3D results [l6j. Frenkel and 
Maguire have exploited the Doi-Edwards theory, valid 
above semi-dilute concentrations, and developed a the- 
ory for the translational motion in 3D 15]. According 



FIG. 10: Fig. 12: Dependence of Dtrs on needles density p. 



to their results C|| scales as D~^^^. Combining this with 
the scaling result Drot ^ and that D^^ is the integral 
of C| I (r) over time they concluded that parallel compo- 
nent of translational diffusion tensor scales as Z?|| ~ p3. 
Moreover, transverse component of diffusion tensor scales 
as D_L ~ p~2 within Doi-Edward theory. Therefore they 
concluded that at high densities D ^ p2 . The MD data 
of Magda et al. [l^ gives the dependence of D± on p 
as D± ~ p"^-^''. Our results in figure (12) is in qualita- 
tive agreement with the result obtained for a 3D gas of 
hard needles 0, El, H. Note that D = ^{D^ + D\\) 
and at high densities it is dominated by the longitudinal 
component. The fitted exponent to the portion p < 5 of 
our result for translational diffusion gives Dtrs ~ p~'^'^^. 
The value 0.93 slightly differs to the predicted exponent 
1 by Enskog theory. This shows dimensionality notably 
affects the system properties. Eventually figure (13) ex- 
hibits the dependence of Drot on p. In 3D based on 
scaling arguments it can be concluded that within the 
Doi-Edward theory Drot scales as p~^. MD simulations 
for a 3D gas of hard needles gives the dependence of Drot 
on density as Drot ~ P"'' with /3 e [1.8,2.2] fl3\. The 
MD results of Magda et al. [11] give the exponent /3 — 1.5 
for 32 < pL^ < 72 and ;3 = 1.89 for the density range 
72 < pL^ < 100 (L is the needle length). We have fitted 
both an algebraic curve Drot ^ p~^ and an exponential 
curve Drot ~ e « to our own data. It turned out that 
P = 3.6 and i = 0.92 with = 0.7828 and 0.9375 for 
the power law and exponential fits correspondingly (ex- 
ponential curve is better fitted). We thus observe a slower 
decrease for Drot in 2D rather than in three. This can 
be explained on the basis of having more pronounced de- 
gree of entanglement and topological constraints among 
needles in two dimensions respect to 3D. 
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FIG. 11: Fig.13: Dependence of Drot on p. 

V. SUMMARY AND CONCLUDING REMARKS 

We have simulated the dynamics of a 2D gas of hard 
needles by event-oriented molecular dynamics. Many of 
the temporal autocorrelation functions both translational 
and angular exhibit a sort of slow dynamics and multi 
step relaxation. The most interesting feature is the exis- 



tence of three regimes in the angular mean squared dis- 
placement. This can be attributed to slow dynamics. 
Our findings show relaxing the translational degrees of 
freedom does not smear out angular slow dynamics. Den- 
sity dependence of translational and rotational diffusion 
coefficients has been obtained and compared to three di- 
mensional results. In 2D dependence of the translational 
diffusion coefficient on density qualitatively resembles to 
3D. Rotational diffusion constant exhibits an algebraic 
decay but with a larger exponent than in 3D. 
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